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ABSTRACT 

To  meet  the  United  States  Government  nuclear  explosion  monitoring  requirements  with  high  confidence,  the  Air 
Force  Technical  Applications  Center  needs  new  and  improved  capabilities  for  analyzing  regional  seismic, 
teleseismic,  and  infrasound  event  data.  Recently,  the  National  Nuclear  Security  Administration  has  decided  to  move 
toward  3D  modeling  to  improve  knowledge  of  the  compressional  and  shear  velocity  structure  and  enable  us  to 
reduce  uncertainty  and  more  accurately  detect,  locate,  and  identify  small  (body  wave  magnitude  mb<4)  seismic 
events.  For  seismically  active  areas,  inaccurate  models  can  be  corrected  using  the  kriging  methodology  and, 
therefore,  it  is  possible  to  detect,  locate,  and  identify  large  events  even  with  limited  resolution  models.  This  is  not 
necessarily  the  case  for  smaller  events,  however,  and  it  is  even  more  of  a  challenge  for  aseismic  regions. 
Furthermore,  interest  on  near-regional  to  local  monitoring  demands  that  we  address  the  Earth’s  heterogeneities  and 
3D  complexities. 

Motivated  by  the  shortcomings  of  existing  single-parameter  inversion  methods  in  accurate  prediction  of  other 
geophysical  parameters,  this  research  was  mainly  focused  on  the  development  and  refinement  of  advanced 
multivariate  inversion  techniques  to  generate  a  realistic,  comprehensive,  and  high-resolution  3D  model  of  the 
seismic  structure  of  the  crust  and  upper  mantle  that  satisfies  multiple  independent  geophysical  datasets.  We  present 
3D  seismic  velocity  models  of  the  crust  and  upper  mantle  beneath  three  different  regions  (northwest  China;  the  East 
Africa  Rift  System;  and  Utah)  resulting  from  the  simulatenous  and  joint  use  of  seismic  body  wave  arrival  times, 
surface  wave  dispersion  measurements,  and  gravity  data.  Our  methodology  represents  a  robust  and  consistent 
compromise  that  fits  the  different  datasets  within  accepted  tolerances.  In  addition  to  obtain  the  optimum  earth  model 
fitting  the  multiple  observations,  we  showed  our  intial  results  towards  an  independent  assessment  of  the  prediction 
capability  of  these  newly  computed  models  using  a  purely  numerical  method  for  wave  propagation  modeling 
(the  Spectral  Element  Method). 
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OBJECTIVES 

The  ultimate  goal  of  this  study  is  to  improve  our  knowledge  of  the  3D  compressional  and  shear  velocity  structure 
and  enable  us  to  reduce  uncertainty  and  more  accurately  detect,  locate,  and  identify  small  (body  wave  magnitude 
mb<4)  seismic  events,  and  therefore  improve  our  capabilities  for  nuclear  explosion  monitoring  (NEM).  This  project 
specifically  improves  seismic  monitoring  technology  through  the  development  and  application  of  advanced 
multivariate  inversion  techniques  to  generate  a  realistic,  comprehensive,  and  high-resolution  3D  model  of  the 
seismic  structure  of  the  crust  and  upper  mantle  that  satisfies  numerous  independent  geophysical  datasets. 


RESEARCH  ACCOMPLISHED 

Inversion  of  surface  wave  dispersion  data  is  a  standard  method  for  determining  3D  shear  velocity  structure  of  the 
crust  and  upper  mantle  of  the  Earth.  On  the  one  hand,  inversion  of  phase  or  group  velocity  dispersion  of  surface 
waves  excited  by  earthquakes  (and  measured  at  relatively  low  frequencies)  has  revealed  shear  wavespeed  variations 
at  wavelengths  upward  of  300  km  (e.g.,  Ritzwoller  and  Levshin,  1998;  Huang  et  al.,  2003;  Lebedev  and  Van  der 
Hilst,  2008;  Yi  et  ah,  2008).  On  the  other  hand,  ambient  noise  tomography  -  with  surface  wave  Green’s  functions 
estimated  from  cross  correlation  of  seismic  ambient  noise  -  has  been  used  to  image  crustal  Vs  variation  with  a 
lateral  resolution  upward  of  100  km  either  on  regional  or  on  sub-continental  scales  (e.g.,  Zheng  et  ah,  2008;  Yang  et 
ah,  2010;  Yao  et  ah,  2010).  Body  wave  travel-time  tomography  and  surface  wave  tomography  each  have  their 
specific  strengths  and  weaknesses.  Travel  time  tomography  can  yield  higher  resolution  in  regions  of  dense  path 
coverage,  and  it  generally  has  excellent  lateral  resolution  beneath  regions  of  high  seismic  activity  or  dense  station 
coverage.  In  regions  of  low  seismicity  and  sparse  station  distribution,  however,  the  shallow  subsurface  cannot  be 
resolved  adequately  by  direct  P  or  S  travel  times.  In  contrast,  surface  wave  data  generally  yield  better  radial 
resolution  and  have  better  potential  for  resolving  shallow  mantle  structure  beneath  regions  that  are  aseismic  or 
which  are  void  of  seismograph  stations.  To  benefit  from  the  complementary  sampling  of  different  seismic  datasets 
(such  as  body  and  surfaces  waves),  multivariate  inversion  techniques  should  be  considered.  This  could  be  achieved 
by  full  waveform  inversion  (e.g.,  Tromp  et  al.,  2005;  Zhao  and  Jordan,  2006),  with  the  implicit  consideration  of 
finite  frequency  kernels  computed  in  heterogeneous  media,  but  the  huge  computational  cost  of  such  an  approach 
makes  it  as  yet  impractical  for  routine  implementation  in  an  operational  setting.  During  this  project  and  with  such 
operational  aspects  in  mind,  we  focused  on  approximations  to  full  wave  propagation  which  are  computationally 
efficient  and  still  sufficiently  accurate  for  the  goal  of  routine  earthquake  location  and  nuclear  monitoring. 

Gravity  measurements  can  provide  constraints  on  spatial  variations  in  (mass)  density  of  rock  in  the  subsurface,  but 
like  any  other  potential  field  method  interpretation  of  gravity  anomalies  is  plagued  by  substantial  ambiguity.  Indeed, 
weak  and  broad  structures  in  the  shallow  subsurface  can  produce  the  same  gravity  signal  (at  the  surface)  as  a  small, 
strong  density  anomaly  at  a  larger  depth.  Using  an  empirical  relationship  between  velocity  and  density,  Maceira  and 
Ammon  (2009)  combined  surface  wave  dispersion  and  gravity  observations  into  a  single  inversion  to  obtain  a  self- 
consistent  high-resolution  3D  shear  velocity  and  density  model  with  increased  resolution  of  shallow  geologic 
structures.  For  a  study  of  Tarim  Basin  (western  China)  they  used  gravity  data  from  the  GRACE  satellite  mission 
(Tapley  et  al.,  2005)  along  with  high-resolution  surface  wave  slowness  tomographic  maps  (Maceira  et  al.,  2005), 
and  the  3D  velocity  model  obtained  from  their  joint  inversion  fits  simultaneously  both  data  sets.  Encouraged  by 
these  results  and  motivated  by  the  shortcomings  of  existing  single-parameter  inversion  methods  in  accurate 
prediction  of  other  geophysical  parameters  (e.g.  Julia  et  al.  2000),  we  have  developed  an  advanced  multivariate 
inversion  technique  to  generate  a  realistic,  comprehensive,  and  high-resolution  3D  model  of  the  seismic  structure  of 
the  crust  and  upper  mantle  that  satisfies  multiple  independent  geophysical  datasets. 

Our  final  algorithm  and  code  JointTomoDD  is  a  modification  of  the  Maceira  and  Ammon  (2009)  joint  inversion 
code,  in  combination  with  the  regional  version  of  the  double-difference  (DD)  tomography  program  tomoDD  (Zhang 
and  Thurber,  2003,  2006),  with  a  fast  LSQR  solver  operating  on  the  gridded  values  jointly.  The  model 
representation  is  a  combination  of  columns  of  rectangular  prisms  (for  gravity  forward  modeling)  embedded  in  a  grid 
whose  nodes  are  located  at  the  center  of  each  prism.  In  a  simplified  manner,  the  system  of  equations  to  be  inverted 
can  be  written  in  the  form: 

WtGtm  =  dt ;  model  includes  compressional  and  shear- wave  velocities 

WgGgm  =  dg ;  model  includes  compressional  and  shear- wave  velocities 
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WsGsm  =  ds ;  model  includes  shear-wave  velocities  only; 

where  wh  wg,  and  are  weighting  parameters  for  seismic  arrival  times,  gravity,  and  surface  wave  dispersion  data 
respectively,  that  equalize  the  three  data  sets.  First-order  smoothing  and  damping  are  also  appplied  to  the  model. 
Please  refer  to  Maceira  et  al.  (2009)  for  further  details. 


East  Africa  Rift  System 

Knowledge  of  lithospheric  structure  is  of  importance  for  understanding  East  Africa's  geodynamic  evolution  and  for 
addressing  broader  questions  about  the  causes  of  continental  breakup.  Though  recent  investigations  have  yielded 
improved  characterizations  of  the  rift  zone,  many  uncertaintites  remain.  For  example,  the  basalt  dominated 
magmatism  in  East  Africa  has  been  explained  by  both  one  deep  mantle  plume  (Ebinger  and  Sleep,  2001)  and  two 
plumes  (Rogers  et  al.,  2000).  Magma-assisted  rifting  in  the  northern  Main  Ethiopian  Rift  contrasts  with  fault- 
controlled  extension  further  south  (Beutel  et  al.,  2010),  but  in  many  cases  the  extent  of  lithospheric  thinning  and 
temporal  and  spatial  evolution  of  rifting  remain  unclear  (e.g.,  Ebinger,  1989).  Key  to  resolving  such  issues  are  better 
constrained  seismic  models.  We  implemented  JointTomoDD  in  this  region.  Benefits  of  our  joint  inversion  approach 
appear  pronounced  when  working  with  regions  of  strong  lateral  contrast  as  found  in  central  Asia  (Maceira  and 
Ammon,  2009).  In  applying  the  joint  inversion  technique  to  East  Africa,  we  solve  for  velocity  structure  in  an  area 
with  less  lateral  heterogeneity  but  great  tectonic  complexity.  To  increase  the  effectiveness  of  the  technique  in  this 
region  we  explore  gravity  filtering  methods  and  test  different  velocity-density  relations  (Modrak  et  al,  2011). 

The  area  for  the  inversion  spans  the 
broad  uplifted  region  from  Ethiopia  at 
one  end  to  Kenya,  Uganda,  and  Tanzania 
at  the  other.  Near  the  northern  boundary 
of  our  study  area,  the  Main  Ethiopia  Rift 
meets  the  incipient  Red  Sea  and  Gulf  of 
Aden  spreading  ridges.  At  the  opposite 
end,  the  rift  system  splits  into  distinct 
western  and  eastern  branches,  which 
largely  sidestep  the  Archean  Tanzania 
%  craton.  Recent  inversions  of  East  Africa 
have  employed  body  waves  (e.g.,  Bastow 
et  al.,  2005;  Benoit  et  al.,  2006),  surface 
waves  (e.g.,  Knox  et  al.,  1998; 
Weeraratne  et  al.,  2003),  receiver 
functions,  or  some  combination  of  these 
(e.g.,  Julia  et  al.,  2005;  Keranen  et  al., 
2009).  Although  useful  comparisons  can 
be  drawn  between  the  Ethiopian  and 
Tanzanian  portions  of  the  rift  system, 
most  tomographic  studies  to  date  have 
focused  exclusively  on  one  section  or  the 
other.  The  current  inversion,  in  contrast, 

is  carried  out  over  a  wider  area  than  most 
previous  studies,  allowing 
straightforward  comparison  between 
these  two  distinct  portions  of  the  rift 
system. 

Fundamental-mode  Rayleigh  wave  group  velocity  estimates  with  periods  from  7  s  to  150  s  were  obtained  from 
Pasyanos  and  Nyblade  (2007)  for  the  inversion.  Though  less  detailed  than  images  from  local  seismic  arrays  (e.g., 
Prodehl  et  al.,  1997),  these  estimates  span  a  broader  spatial  and  period  range.  Gravity  data  for  the  inversion  were 
derived  from  the  Gravity  Recovery  and  Climate  Experiment  (GRACE)  (Tapley  et  al.,  2005).  We  implemented  a 
method  to  increase  the  usefulness  of  gravity  data  by  filtering  the  Bouguer  anomaly  map.  Though  commonly  applied 


40  km 

30D  35D  (CO 


60  km 

30D  35D  fCD 


80  km 

30  □  350  ®m 


□§□10  » 


□4  OB  0  ■ 


□  10  1 


Figure  1.  Horizontal  8Vs/Vs  cross-sections  at  various  depths 

through  the  3D  velocity  model  obtained  from  the  joint 
inversion.  Percent  values  are  with  respect  to  mean 
velocities  shown  in  the  corner  of  each  cross-section. 
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in  gravity  forward  modeling  (e.g.,  Simiyu  and  Keller,  1997),  such  techniques  have  not  to  our  knowledge  been  used 
in  previous  joint  inversion  studies  (e.g.,  Lees  and  VanDecar,  1991;  Zeyen  and  Achauer  1997;  Tiberi  et  al.,  2003). 
Different  tests  suggest  that  addition  of  unfiltered  gravity  data  contributes  little  in  the  way  of  distinguishing  between 
features  at  different  depths.  Rather  than  improving  resolution  at  shallow  depths  as  desired,  features  in  the  unfiltered 
gravity  data  are  smeared  into  the  mantle.  Although  filtering  removes  potentially  useful  information  on  mantle 
structure,  the  remaining  short-wavelength  signal  can  be  assigned  with  greater  reliability  within  the  crust,  avoiding 
the  mutually  degrading  effects  of  smearing  between  crust  and  mantle.  To  remove  the  long-wavelength  components 
from  the  Bouguer  gravity  map  we  follow  Tessema  and  Antoine  (2004),  who  use  an  upward  continuation  method 
and  demonstrate  correlation  with  crustal  geology. 


Figure  1  shows  the  3D  S-wave  velocity 
model  obtained  from  the  joint  inversion. 

The  low- velocity  anomaly  beneath  Ethiopia 
is  among  the  most  prominent  features.  The 
anomaly  appears  most  conspicuous  at  -60 
km  depth  beneath  Afar  and  continues 
southward  along  the  Main  Ethiopian  Rift  at 
greater  depths.  Although  low  velocities 
beneath  Ethiopia  appear  pronounced  up  to 
-140  km,  Fig.  1  suggests  the  magnitude  of 
the  anomaly  becomes  somewhat  diminished 
by  -150  km.  Such  a  result  appears  in 
agreement  with  a  number  of  regional 
surface  wave  studies.  For  example,  while 
Ritsema  and  Van  Heijst  (2000)  resolve  a 
low- velocity  anomaly  beneath  Ethiopia 
extending  to  at  least  250  km,  a  decrease  in 
the  magnitude  of  this  anomaly  becomes 
visible  at  around  1 60  km.  Similar 
magnitude  decreases  are  visible  in  the 
results  of  Knox  et  al.  (1998)  and  Pasyanos 
andNyblade  (2007).  Besides  the  low- 
velocity  anomaly  below  Ethiopia, 
prominent  velocity  excursions  also  occur 
below  Tanzania.  In  the  40  km  depth  slice 
we  resolve  lower  shear  velocities  beneath 
the  Tanzania  craton  than  in  the  adjoining 
rift  branches;  at  this  depth,  lower  velocities 
beneath  the  craton  are  readily  explained  by 
the  contrast  between  thick  crust  in  the 
craton  and  shallow  mantle  in  the 
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Figure  2.  Fit-to-data  from  inversion  of  group  velocities  only  and 
from  inversion  of  group  velocities  and  gravity,  (a)  Top: 
Group  velocities  from  a  representative  cell  in  the 
model.  Bottom:  Filtered  Bouguer  anomalies,  (b)  Top: 
Group  velocity  fit  obtained  from  inversion  of  group 
velocities  only.  Bottom:  Gravity  fit  obtained  from 
inversion  of  group  velocities  only,  (c)  Top:  Group 
velocity  fit  obtained  from  joint  inversion.  Bottom: 
Gravity  fit  obtained  from  joint  inversion. 


surrounding  rift  branches.  Beginning  at  50  km,  velocities  in  the  craton  revert  to  higher  values  than  in  the  adjacent 
non-cratonic  terranes;  these  higher  values  persist  to  -120  km.  Finally,  at  -140  km,  a  pronounced  low- velocity 
anomaly  emerges  beneath  the  craton.  This  juxtaposition  of  high  and  low  shear- wave  speeds  between  120-140  km 
depth  appears  consistent  with  the  hypothesis,  discussed  in  detail  by  Weeraratne  et  al.  (2003)  and  Nyblade  et  al. 
(2000),  of  a  hot,  upwelling  plume  impeded  by  cool,  overyling  material  of  the  craton.  Additionally,  our  results  allow 
comparison  between  rift  structure  of  Ethiopia  and  Tanzania.  In  obtaining  data  from  Pasyanos  and  Nyblade  (2007), 
we  use  group  velocities  derived  not  only  from  local  stations  and  events,  but  also  from  stations  and  events  distributed 
across  surrounding  tectonic  plates.  Though  the  resulting  continent-scale  maps  possess  less  detail  than  local-scale 
group  velocity  maps,  their  wider  spatial  coverage  allows  straight-forward  comparison  between  distinct  portions  of 
the  rift  system.  As  a  result,  we  find  that  uppermost  mantle  shear  velocities  beneath  Ethiopia  appear  much  slower 
than  those  beneath  Tanzania.  While  the  presence  of  shallow  low  velocities  beneath  Ethiopia  suggests  that  mantle 
lithosphere  there  has  been  largely  replaced  by  asthenosphere  (e.g.,  Beutel  et  al.,  2010),  the  absence  of  shallow  low 
velocities  beneath  the  southern  rift  branches  is  consistent  with  fault-controlled  extension  in  that  part  of  the  rift 
system.  Finally,  though  a  common  origin  at  greater  depths  is  not  ruled  out,  no  evidence  is  observed  that  the  various 
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low- velocity  anomalies  in  Fig.  1  merge  continuously  above  175  km.  This  possibly  explains  why  igneous  rocks  from 
the  two  low  velocity  zones  are  compositionally  different  (e.g.,  Rogers  et  al.,  2000). 

In  the  inversion  carried  out  in  central  Asia  by  Maceira  and  Ammon  (2009),  addition  of  gravity  data  dramatically 
improved  the  fit  to  the  Bouguer  anomalies  without  significantly  degrading  the  fit  to  the  group  velocities.  Figure  2 
demonstrates  this  result  for  the  current  study  area  as  well.  Although  it  is  well  known  that  problems  of  non¬ 
uniqueness  make  gravity  data  easier  to  match  than  seismic  data,  several  observations  provide  confidence  in  our 
methodology's  robustness.  These  include  the  simultaneous  fit  to  both  data  sets  shown  in  Fig.  2  as  well  as  qualitative 
changes  resulting  from  the  addition  of  gravity  data.  For  example,  compared  with  results  from  the  inversion  of  group 
velocities  only,  the  joint  inversion  methodology  provides  increased  effectiveness  capturing  Moho  depth  at  the 
continental  margin  and  sharper  delineation  of  the  Tanzania  craton.  The  resolved  extent  of  the  high-velocity  cratonic 
region  accords  well  with  previous  tomographic  images  (e.g.,  Fig.  1 1  of  Weeraratne  et  al.,  2003)  as  well  as 
geodynamic  models  suggesting  strain  localization  in  zones  of  weakness  surrounding  the  craton  (e.g.,  Nyblade  and 
Brazier,  2002). 

Utah  Geothermal  Region 

The  Cove  Fort-Sulphurdale  geothermal  area  is  located  in  the  transition 
zone  between  the  Basin  and  Range  to  the  west  and  the  Colorado  Plateau  to 
the  east  (Figure  3).  We  have  collected  various  geophysical  data  around  the 
geothermal  field,  including  gravity  anomalies  (Pan-American  Center  for 
Earth  and  Environmental  Studies  (PACES)  available  at 
http://gis.utep.edu),  seismic  surface  wave  phase  and  group  velocity  maps 
(Yang  et  al.,  2008),  and  seismic  body  wave  arrival  times  that  were 
assembled  from  seismic  waveforms  recorded  by  the  University  of  Utah 
Seismograph  Stations  (UUSS)  regional  network  for  the  past  7  years  and 
the  recent  EArthscope/US Array  phase  data.  Various  geophysical  data  sets 
indicate  that  beneath  the  Cove  Fort-Sulphurdale  geothermal  resource  there 
is  a  strong  anomaly  of  low  seismic  velocity,  low  gravity  and  high  electrical 
conductivity  that  correlates  with  the  high  surface  heat  flow.  This  suggests 
that  there  is  a  heat  source  in  the  crust  beneath  the  geothermal  field.  We 
collected  first  P  arrival  data  from  more  than  6500  earthquakes  in  the  Utah 
region.  Each  event  has  at  least  6  arrivals  for  reliably  determining  its 
location.  We  applied  the  double-difference  seismic  tomography  method 
(Zhang  and  Thurber,  2003)  to  simultaneously  determine  an  initial  velocity 
Figure  3.  Simplfied  tectonic  map  structure  and  earthquake  locations.  On  the  preliminary  regional  seismic 
showing  tectonic  provinces  velocity  map  computed  this  way,  we  can  also  identify  some  other  low 
around  Cove  Fort  velocity  anomalies,  indicating  other  potential  geothermal  prospects.  We 

geothermal  field.  then  applied  our  simultaneous  joint  inversion  methodology  to  produce  a 

better  constrained  velocity  structure  of  the  Utah  area  which  will  be  very 
helpful  for  characterizing  and  exploring  existing  and  potential  geothermal  reservoirs  in  the  area. 


Relative  weighting  between 
gravity  and  surface  waves 


Surface  wave  data  residual 


Relative  weighting  between  Relative  weighting  between 

surface  waves  and  travel  times  gravity  and  travel  times 


Figure  4.  Trade-off  curves  between  datasets  pairs  are  used  to  decide  on  relative  weighting. 
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Joint  inversion  of  seismic  travel  times,  surface  wave  dispersion,  and  gravity  data  represents  a  multiple-objective 
optimization  problem.  Because  it  is  unlikely  that  the  different  “objectives”  (data  types)  would  be  optimized  by  the 
same  parameter  choices,  some  trade-off  between  the  objectives  is  needed.  Figure  4  shows  an  example  of  finding  the 
relative  weightings  between  the  multiple  data  types  through  a  trade-off  analysis  of  data  residuals.  As  a  result,  the 
final  model  will  optimally  fit  the  different  datasets. 


Travel  Time  RMS  residual:  1 47  ms  Travel  Time  RMS  residual:  1 89  ms 

Gravity  RMS  residual:  20  mgal  Gravity  RMS  residual:  0.23  mgal 


Figure  5.  Compressional-wave  velocity  model  at  constant  depth  slices  using  (a)  seismic  travel  times 
alone  and  (b)  joint  inversion  of  body  waves  travel  times  and  gravity.  (Velocity  units:  km/s). 

Figures  5  and  6  show  different  depth  slices  through  the  computed  model  for  compressional  and  shear- wave  velocity 
respectively.  Figure  5  shows  the  comparison  of  the  Vp  model  obtained  from  travel  time  arrivals  only  with  that 
obtained  using  seismic  arrivals  together  with  gravity  anomaly  information.  Figure  6  shows  the  comparison  of  the  Vs 
model  obtained  from  the  joint  inversion  of  surface  wave  and  gravity  data  with  that  obtained  using  all  three  data 
types.  The  latter  model  fits  the  three  data  sets  well  and  shows  better  definition  of  the  velocity  anomaly  associated 
with  the  transition  from  the  Basin  and  Range  to  the  west  to  the  Colorado  Plateau  to  the  east. 


Figure  6.  Shear-wave  velocity  model  at  a  depth  of  17  km  obtained  from  (A)  surface  wave  data,  (B)  surface 
and  gravity  data,  (C)  surface  and  travel  time  data,  and  (D)  surface,  gravity,  and  travel  time  data. 
(Velocity  units:  km/s). 
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Northwest  China 


Our  final  goal  is  to  generate  a  3D  realistic  and  comprehensive  model  of  the  crustal  and  upper  mantle  seismic 
structure  underneath  northwest  China,  an  area  of  prime  importance  for  the  nuclear  explosion  monitoring  program. 
We  have  obtained  a  model  from  the  joint  inversion  of  surface  waves  dispersion  measurements  (Maceira  et  ah, 
2005),  teleseismic  P-wave  receiver  functions  (Ammon  et  al.,  2004),  and  gravity  anomalies  (Tapley  et  al.,  2005). 
This  model  fits  simultaneously  all  the  data  sets  offering  a  compromise  between  fitting  the  three  data  sets  (Figure  7). 

We  are  now  in  the  process  of  refining  a  new  3D  model  obtained  by  incorparating  a  fourth  dataset  in  the  inversion 
process.  Body  waves  (P  and  S)  travel  times  are  gathered  from  the  LANL  Knowledge  Database.  Figure  8  shows 
preliminary  results  which  are  in  good 

agreement  with  geological  and  tectonic  observations 

knowledge  of  the  area.  Final  results  will  be 
shown  during  the  Monitoring  Research 
Review  in  September. 


RESULTS  FROM  INVERSION 
OF  SURFACE  WAVES  AND 
RECEIVER  FUNCTIONS 


RESULTS  FROM  JOINT 
INVERSION  OF  THREE 
DATA  SETS 


Figure  7.  Fit  to  the  different  datasets  -  P- 
wave  receiver  functions  (top  row),  surface 
wave  dispersion  (middle  row),  gravity 
(bottom  row)  -  from  the  inversion  of 
surface  waves  and  receiver  functions 
(middle  column)  and  from  the  joint 
inversion  of  the  three  datasets  (right 
column). 
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Figure  8.  Vs  model  depth  slices  from  joint  inversion  of  body  waves,  surface  waves,  and  gravity  anomalies. 
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Model  Validation 

To  address  the  need  of  near-regional  to  local 
nuclear  explosion  monitoring,  several  research 
institutions  have  been  focused  on  inferring  the 
best  resolution  possible  images  of  the 
underground  solid  Earth.  For  the  last  three  years, 
LANL  has  been  developing  and  applying 
advanced  multivariate  inversion  techniques  to 
generate  realistic  and  high-resolution  3D  models 
of  the  seismic  structure  that  satisfies  numerous 
independent  geophysical  datasets.  Our  more 
complete  models  obtained  by  simultaneous  joint 
inversion  of  disparate  data  sets  produce  better 

predictions  and  minimize  the  differences  between 
observations  and  predictions.  However,  while  our 
best  current  inversion  techniques  are  providing 
3D  velocity  models  with  the  best  resolution  ever, 
they  don’t  provide  any  absolute  assessment  of  the 
model  uncertainty.  Theoretically,  this  will  require 
testing  the  entire  range  of  possible  values  for  each  parameter  to  get  a  complete  “picture”  of  the  solution  space;  in 
addition  to  questioning  the  fidelity  of  the  imaging  method  (in  this  case,  the  seismic  wave  propagation  scheme  used 
to  predict  waveforms). 

Geophysical  model  validation  is  typically  done  using  resolution  tests  that  assume  the  imaging  theory  used  is 
accurate  and  thus  only  considers  the  impact  of  the  data  coverage  for  resolution.  We  are  taking  a  more  rigorous 
approach  to  model  validation  via  full-waveform  propagation.  LANL  High  Performance  Computing  resources  allow 
us  to  perform  accurate  3D  modeling  of  wave  propagation  through  these  models  using  the  Spectral  Element  Method 
(SEM).  This  recently  used  in  seismology  technique  makes  no  assumptions  about  the  theory  used  to  generate  the 
models  but  require  substantial  computational  resources.  SEM  is  a  particular  case  of  continuous  Galerkin  method 
with  optimized  efficiency  because  of  its  use  of  high  order  and  tensorized  basis  functions  (Komatitsch  and  Tromp, 
1999).  It  can  handle  very  distorted  elements  (Oliveira  &  Seriani,  2011)  imposed  by  complex  geophysical  models. 
The  parallel  implementation  of  SEM  utilizes  domain  partition  to  partition  the  elements  amongst  the  compute  nodes. 
Current  implementation  incorporates  3D  topography  of  seismic  interfaces,  anisotropy  and  attenuation. 

We  are  currently  and  systematically  computing  the  misfit  between  predicted  and  actually  observed  waveforms  for 
the  different  Earth  models  generated  in  the  project.  Due  to  operational  delays  arising  from  the  emergency  LANL 
closure  and  incorrect  shoutdown  of  the  HPC  resources,  we  will  present  these  results  during  the  coming  MRR. 

Figure  9  shows  the  idea  of  this  model  validation  but  with  the  DNA09-Berkeley  model.  The  tests  performed  for  this 
model  and  a  moderate-sized  event  on  the  Pacific  Northwest  show  no  perceptible  difference  between  models 
obtained  with  two  different  imaging  techniques  (finite-frequency  ray-theoretical)  at  intermediate  periods. 
Differences  star  to  appear,  however,  at  higher  frequencies. 

CONCLUSIONS  AND  RECOMMENDATIONS 

This  is  the  last  year  of  a  three-year  project  to  map  the  three-dimensional  (3D)  seismic  structure  of  the  crust  and 
upper  mantle  using  simultaneous  joint  inversion  of  surface  wave  dispersion,  gravity,  receiver  function,  and  travel 
time  observations.  We  have  successfully  accomplished  our  goal  of  developing  an  algorithm  and  corresponding 
computer  codes  for  advanced  multivariate  inversion  for  Earth  structure.  We  have  dealt  with  multiple  challenges  of 
multivariate  approaches  such  as  relationship  between  independent  variables  in  the  inversion  scheme  and  relative 
weighting  of  disparate  datasets.  We  have  also  learned  that  besides  enhancing  resolution  at  short  wavelengths,  use  of 
filtered  gravity  anomalies  may  help  distinguish  between  anomalies  at  different  depths.  We  are  now  refining  and 
validating  our  3D  models  for  different  areas  around  the  globe.  Our  more  rigorous  approach  to  model  validation  via 
full-waveform  propagation  will,  in  the  near- future,  allow  us  to  quantify  model  uncertainties  and  map  them  into  the 
uncertainty  in  seismic  source  parameters. 


Full  bandwidth  (>  10s)  50-200  s 
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Correlation  coefficient  Correlation  coefficient 

Figure  9.  Cross  correlation  coefficient  between  synthetics 
computed  from  DNA  model  with  finite-frequency 
and  ray-theoretical  approaches;  (left)  considering 
high  frequencies,  (right)  only  periods  50-200  s. 
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